Basis set based orbital generation with NWChem

Basis set based orbital generation with NWChem#

Parameters#

distance = 2.5
iteration_energies = []
iterations = 6
molecule_name = "beh2"
box_size = 50.0
wavelet_order = 7
madness_thresh = 0.0001
basisset = 'cc-pvdz'

Run NWChem calculation#

Create NWChem input and run NWChem calculation. If the MadPy devcontainer or the singularity image is used, NWChem is already installed. Otherwise, NWChem has to be installed and the path has to be adjusted.

import subprocess as sp
nwchem_input = '''
title "molecule"
memory stack 1500 mb heap 100 mb global 1400 mb
charge 0  
geometry noautosym nocenter
  Be 0.0 0.0 0.0
  H 0.0 0.0 ''' + distance.__str__() + '''
  H 0.0 0.0 ''' + (-distance).__str__() + '''
end
basis  
  * library ''' + basisset + '''
end
scf  
 maxiter 200
end   
task scf  
'''

with open("nwchem", "w") as f:
    f.write(nwchem_input)

programm = sp.call("/opt/conda/bin/nwchem nwchem", stdout=open('nwchem.out', 'w'), stderr=open('nwchem_err.log', 'w'), shell = True)

Convert NWChem AOs and MOs to MRA-Orbitals#

Read the atomic orbitals (AOs) and molecular orbitals (MOs) from a NWChem calculation and translate them into multiwavelets.

import madpy as mad
red = mad.RedirectOutput("NWChem_to_MRA.log")
converter = mad.NWChem_Converter(box_size, wavelet_order, madness_thresh)

converter.Read_NWChem_File("nwchem")
aos = converter.GetNormalizedAOs()
mos = converter.GetMOs()

del converter
del red

Plotting of orbitals#

First, a molecule object with the same geometry as in the NWChem calculation is created and then orbital(s) (here MO 5) are exported to a .cube file.

molecule = mad.molecule()
molecule.add_atom(0, 0, 0, "Be")
molecule.add_atom(0, 0, distance, "H")
molecule.add_atom(0, 0, -distance, "H")

red = mad.RedirectOutput("plot.log")
plotter = mad.Plot(box_size, wavelet_order, madness_thresh)
plotter.cube_plot("orbital", mos[5], molecule, datapoints=101)
del plotter
del red

The .cube file can then be loaded and plotted with py3Dmol, for example.

import py3Dmol
orbdata = open("orbital.cube", "r").read()
v = py3Dmol.view()
v.addVolumetricData(orbdata, "cube", {'isoval': -0.001, 'color': "red", 'opacity': 0.75})
v.addVolumetricData(orbdata, "cube", {'isoval': 0.001, 'color': "blue", 'opacity': 0.75})
v.addModel(orbdata, "cube")
v.setStyle({'sphere':{}})
v.zoomTo()
v.show()

3Dmol.js failed to load for some reason. Please check your browser console for error messages.

Contents